set.seed(123)

c <- 1
n <- 100000

f_exp <- function(x, y)
  exp(-x * y)

M <- 1.1

samples <- matrix(0, nrow = 0, ncol = 2)

while (nrow(samples) < n) {
  x_prop <- runif(1000, 0, c)
  y_prop <- runif(1000, 0, c)
  z <- runif(1000, 0, M)
  
  
  accept <- z < f_exp(x_prop, y_prop)
  
  accept_point <- cbind(x_prop[accept], y_prop[accept])
  
  samples <- rbind(samples, accept_point)
}

samples <- samples[1:n, ]

colnames(samples) <- c('x', 'y')


plot(
  samples[, 'x'],
  samples[, 'y'],
  xlab = 'x',
  ylab = 'y',
  main = "Scatter Plot of Samples from p(x,y) Θ e^(-xy)",
  pch = 19,
  col = rgb(0, 0, 1, alpha = 0.5),
  xlim = c(0, c),
  ylim = c(0, c)
)


gibbs_smapler <- function(n_iter, x0, y0) {
  samples <- matrix(0, nrow = n_iter, ncol = 2)
  
  colnames(samples) <- c('x', 'y')
  
  x <- x0
  y <- y0
  
  
  for (i in 1:n_iter) {
    x <- rexp(1, rate=y)
    y <- rexp(1, rate=x)
    
    samples[i,] <- c(x,y)
  }
  
  return(samples)
}


n_iter <- 50000
x0 <- 1
y0 <- 1

samples <- gibbs_smapler(n_iter, x0, y0)

plot(samples[,1], samples[,2],
     xlab='x', ylab='y',
     main='Giibs Sampling for p(x,y) Θ exp(-xy)',
     pch=20, col=rgb(0,0,1, alpha=0.1)
     )

rtrunc_exp <- function(n, rate, c){
  u <- runif(n)
  
  return(-log(1- u* (1-exp(-rate * c)))/rate)
}

gibbs_sampler <- function(n_iter, c, x0, y0, burn_in = 0){
  samples <- matrix(0, nrow=n_iter, ncol=2)
  colnames(samples) <- c('x', 'y')
  
  x <- x0
  y <- y0
  
  for(xi in 1 : n_iter){
    x <- rtrunc_exp(1, rate=y, c=c)
    y<- rtrunc_exp(1, rate=x, c=c)
    
    samples[xi, ] <- c(x,y)
  }
  
  if(burn_in > 0){
    samples <- samples[-(1:burn_in),]
  }
  
  return(samples)
}


n_iter <- 100000
c <- 5
x0 <- 1
y0 <- 1

burn_in <- 1000


samples <- gibbs_sampler(n_iter, c, x0, y0, burn_in)


plot(samples[,1], samples[,2],
     xlab='x', ylab='y',
     main=paste("Gibbs samples for p(x,y) = exp(-xy), 0<x,y<", c),
     pch=20, col=rgb(0,0,1,alpha=0.1),
     xlim=c(0,c),ylim=c(0,c))

par(mfrow=c(1,2))

hist(samples[,1], breaks=50, main='Marginal Distribution of x',
     xlab='x', col='lightblue', probability = TRUE, xlim=c(0,c))
hist(samples[,2], breaks=50, main='Marginal Distribution of y',
    xlab='y', col='lightgreen', probability = TRUE, xlim=c(0,c))


n_size <- 50
x <- sample(1:100, n_size)
y <- sample(1:100, n_size)


plot_data <- samples[1:n_size,]

colnames(plot_data) <- c('x', 'y')



n_plot <- length(x)
indices <- 1:n_plot

plot(plot_data[,1], plot_data[,2], type='n', xlab='x', ylab='y', main='Movement Path of Points')

points(plot_data[,1], plot_data[,2], pch=20, col='blue')

text(plot_data[,1], plot_data[,2], labels = indices, pos=3, cex=0., col='black')

# lines(plot_data[,1], plot_data[,2], col='blue', lty=2)

for(i in 1 : (n_plot-1)){
  x_start <- plot_data[i,1]
  y_start <- plot_data[i,2]
  
  x_end <- plot_data[i+1,1]
  y_end <- plot_data[i+1,2]

  segments(x0=x_start, y0=y_start, x1=x_end, y1=y_start, col='blue', lty=2)
  segments(x0=x_end, y0=y_start, x1=x_end, y1=y_end, col='blue', lty=2)

}

# for(i in 1:(n_plot-1)){
#   arrows(plot_data[i,1], plot_data[i,2], plot_data[i+1,1], plot_data[i+1,2],
#          col="blue", length=0.1, lwd=1)
# }











